%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                          Counterfactual Model                           %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [g, LFN, AvFirmSize, CostMarkup_A, CostMarkup, EntryRate, Lam, M, prod_share, LPNwMLam, k_func_lam, k_func_mu, ExitRate, xProds, ProductDist, Firms_10y, FirmAge_dens, Exit_age, Esales] = Counterfactual(p)

% ---------------------------- Description --------------------------------

    % This function computes the main objects of the model given a vector of
    % parameters.

%--------------------------------------------------------------------------
%                            Preliminaries                                %
%--------------------------------------------------------------------------

    fprintf('------------------------- Counterfactual stage started ----------------------\n')

    % ---------------------------------------------------------------------
    %                 1. Solving for Direct Parameters                    %
    % ---------------------------------------------------------------------

    % Rate of incumbent product creation - see Proposition (M-3)
    p.xstar = (p.phix/(p.zeta*p.phie))^(1/(p.zeta-1));

     % Creative destruction rate - see Proposition (M-3)
    p.tau = p.alpha*(p.eta+p.d0)/(1-p.alpha);

    % Flow entry rate - see Proposition (M-3)
    p.z = 1/p.alpha*p.tau-p.xstar;

    % Net rate of product accumulation - see Section (RP-2.2.2)
    p.psi = p.xstar - (p.tau+p.d0);

    % Average Efficiency gain of new products raised to sigma-1 - see Equation (M-5)
    p.qbar = p.alpha*p.l^(p.s-1)+(1-p.alpha)*p.wbar;

    % The CES markup level - see Section (M-2)
    p.mu = p.s/(p.s-1);

    %------------------------- Computing Moments --------------------------

    % Compute the Aggregate Productivity Growth rate - see Equation (RP-1)
    g = p.eta/(p.s-1) + p.In + p.tau/(p.alpha*(p.s-1))*(p.qbar-1);

    % Computes the numerator of Firm Age Density - see Equation (RP-4)
    Firms_function = @(a) p.psi*exp((p.psi-p.eta).*a)./(p.psi - p.xstar*(1-exp(p.psi.*a)));

    % Integrates the numerator over firm age for denominator - see Equation (RP-4)
    Firms = integral(Firms_function, 0, Inf);

    % Computes firm age density as a function - see Equation (RP-4)
    FirmAge_dens = @(a) Firms_function(a) / Firms;

    % Computes the average number of products per firm - see Equation (RP-5)
    xProds = (p.z * Firms)^(-1);

    % Computes the aggregate entry rate - see Equation (RP-6)
    EntryRate = p.z * xProds;

    %---------------------- Product Distribution --------------------------

    % Preallocate a vector for the unconditional product distribution
    ProductDist = zeros(1, p.maxprods);

    % Compute the Gamma function - see Equation (RP-3)
    GammaFunc = @(a) p.xstar*(1-exp(p.psi*a))./(p.xstar*(1-exp(p.psi*a))-p.psi);

    % Computes the unconditional product probabilities recursively...
    for n = 1:p.maxprods

        % ... by defining the integrand - see Equations (RP-2) and (RP-7)
        Integrand = @(a) (1-GammaFunc(a)) .* GammaFunc(a).^(n-1) .* FirmAge_dens(a);

        % ... and summing over firm age - see Equations (RP-7)
        ProductDist(n) = integral(Integrand, 0, Inf);

    end

    % Clear the integrand function
    clear Integrand

    % Compute the Exit Rate implied by the model - see Equation (RP-52)
    ExitRate = (p.tau + p.d0) * ProductDist(1);

    % Share of firms older than 10 years - see Equation (RP-56)
    Firms_10y = integral(FirmAge_dens, 10, Inf);

    % ---------------------------------------------------------------------
    %                  2. Solving for Average Firm Size                   %
    % ---------------------------------------------------------------------  

    % For an overview of this part of the code, please refer to Section 2.3
    % of the Replication Package

    %----------------------------------------------------------------------

    % The constant C shows up in the differential equation for k(Delta) - see Section (RP-2.3.1)
    ConstC = p.rho + p.d0 + (p.eta+p.d0)*(p.qbar/(1-p.alpha)-1);

    % Boundary condition delivers k(s/s-1) - see Section (RP-2.3.1)
    k_func_mu = 1/(ConstC*(p.s-1))*(p.s/(p.s-1))^(-p.s);

    % Constant of integration for the diff. equation - see Equation (RP-12)
    ODE_const = (1/(ConstC*(p.s-1))-p.s/(ConstC*(p.s-1)+p.In*(p.s-1)^2)+1/(ConstC+p.In*p.s))*(p.s/(p.s-1))^(-p.s-ConstC/p.In);

    % Solution of differential equation to obtain k(lambda) - see Equation (RP-11)
    k_func_lam = p.l^(-p.s)*(p.l/(ConstC+p.In*(p.s-1))-1/(ConstC+p.In*p.s)) + ODE_const*p.l^(ConstC/p.In);
    
    % -------------- Joint Distribution of Gaps and Productivity ----------

    % Scale parameter for Pareto distribution - see Equation (RP-16)
    p.wmin = (p.wbar*(1-(p.s-1)/p.varrho))^(1/(p.s-1)); 

    % Growth rate of average efficiency - see Section (RP-2.3.2)
    p.gQ = (p.tau/p.alpha)*(p.qbar-1)/(p.s-1) + p.In;

    % Computes joint distributions and misallocation terms - see Section (RP-2.3.2) 
    [FNC, F1, M, Lam, Deltagrid, p.qgrid_adj] = Joint_Dist(p);

    % Compute ls^P (times Misallocation terms) - see Equation (RP-9)
    LPNwMLam = (1 - (p.zeta-1)*p.xstar/(p.zeta*(p.rho+p.tau+p.d0)))/(p.phie*(p.alpha*k_func_lam*p.l^(p.s-1)+(1-p.alpha)*k_func_mu*p.wbar));

    % Compute ls^P by removing the misallocation terms accordingly
    LPN = LPNwMLam * M^(p.s-1)*Lam^(p.s);

    % With ls^P, we can compute l = L/N - see Equation (RP-8)
    LFN = LPN + ((p.eta+p.d0)/(1-p.alpha) - (p.zeta-1)*p.xstar/p.zeta)/p.phie;

    % Having ls^P and l, compute production share s^P on the BGP
    prod_share = LPN/LFN;

    % Compute the Average Firm Size - see Equation (RP-22)
    AvFirmSize = LFN * xProds;

    % ---------------------------------------------------------------------
    %                 3. Mark-up and Sales growth                         %
    % ---------------------------------------------------------------------

    % For an overview of this part of the code, please refer to Section 2.4
    % of the Replication Package

    %----------------------------------------------------------------------

    % Conditional distribution of product age given firm age and products - see Section (RP-2.4.4)
    [AgeDist, avec, Pmatrix, p0] = Conditional_Dist(p);

    % Take differences to obtain the associated density
    TimePeriods = p.maxTime/p.periodlength + 1;
    AgeDensity = zeros(TimePeriods, TimePeriods, p.maxprods);
    AgeDensity(2:end,:,:) = diff(AgeDist, 1, 1);
    AgeDensity(1,:,:) = AgeDist(1,:,:);

    % Conditions conditional product probabilities on firm survival
    CondProbs =  Pmatrix ./ repmat((1-p0), p.maxprods, 1);

    % ------------------------ a. Mark-ups --------------------------------

    % Markup by product age for varieties with competitors - see Equation (RP-24)
    markupbyage_C = min(p.l * exp(p.In*avec),(p.s/(p.s-1)));

    % Expected markup by product age - see Equation (RP-25)
    markupbyage = p.alpha*markupbyage_C + (1-p.alpha)*p.mu;

    % Expected markup given firm age and number of products - see Equation (RP-23)
    Emarkup_An = permute(sum(AgeDensity .* repmat(markupbyage', 1, TimePeriods, p.maxprods), 1),[3 2 1]);
    
    % Expected markup given firm age - see Equation (RP-23)
    Emarkup = sum(Emarkup_An .* CondProbs, 1);

    % Finds the index on age grid associated with a firm of 10 years
    [~, ind10] = min(abs(avec-10));

    % Computes markup growth by age 10 - see Equation (RP-30)
    MarkupGrowth = Emarkup(ind10) - Emarkup(1);

    % ------------------------ b. Sales -----------------------------------

    % Expected product-level sales by product age - see Equation (RP-26)
    salesbyage = exp((p.s-1)*(p.In-p.gQ).*avec) / (M^(p.s-1)*Lam^(p.s-1)) .* (p.alpha*p.l^(p.s-1)*markupbyage_C.^(1-p.s) + (1-p.alpha)*p.wbar*p.mu^(1-p.s));

    % Expected product-level sales given firm age and number of products - see Equation (RP-23)
    Sales_An = permute(sum(AgeDensity .* repmat(salesbyage', 1, TimePeriods, p.maxprods), 1),[3 2 1]);

    % Vector with number of products to turn product level into firm level
    prodmat = repmat((1:p.maxprods)',1,TimePeriods);

    % Expected firm-level sales given firm age - see Equation (RP-23)
    Esales = sum(Sales_An .* CondProbs .* prodmat, 1);

    % Computes sales growth as a function of firm age
    Esales = log(Esales) - log(Esales(1));

    % Computes firm-level sales growth by age 10 - see Equation (RP-30)
    SalesGrowth = Esales(ind10);

    % ------------------ c. Cost Weighted Markup --------------------------

    % Expected cost-weighted markup by product age - see Equation (RP-27)
    costmarkupsbyage = (p.alpha*p.l^(p.s-1)*markupbyage_C.^(1-p.s) + (1-p.alpha)*p.wbar*p.mu^(1-p.s)) ./ (p.alpha*p.l^(p.s-1)*markupbyage_C.^(-p.s) + (1-p.alpha)*p.wbar*p.mu^(-p.s));

    % Expected cost-weighted markup given firm age and number of products - see Equation (RP-23)
    CostMarkup_An = permute(sum(AgeDensity .* repmat(costmarkupsbyage', 1, TimePeriods, p.maxprods), 1),[3 2 1]);

    % Expected cost-weighted markup given firm age - see Equation (RP-23)
    CostMarkup_A = sum(CostMarkup_An .* CondProbs, 1);

    % Integrate over firm age to get aggregate markup - see Equation (RP-30)
    Integrand = @(a) interp1(avec, CostMarkup_A, a, "linear", CostMarkup_A(end)) .* FirmAge_dens(a);
    CostMarkup = integral(Integrand, 0, Inf);

    % Clear the remaining variables
    clear markupbyage_C markupbyage Emarkup_An ind10 salesbyage Sales_An prodmat costbyage markupcostbyage costmarkupsbyage CostMarkup_An Integrand

    % ---------------------------------------------------------------------
    %                   4. Employment Distribution                        %
    % ---------------------------------------------------------------------

    % For an overview of this part of the code, please refer to Section 2.5
    % of the Replication Package

    %----------------------------------------------------------------------

    % Create a matrix with Firm Age density with expanded dimensions
    FirmAge_matrix = repmat(FirmAge_dens(avec), TimePeriods, 1, p.maxprods);

    % Obtain conditional distribution of product age given number of products - see Section (RP-2.5)
    AgeDist_N = permute(trapz(avec, AgeDist .* FirmAge_matrix, 2), [1, 3, 2]);

    % Take differences from distribution to get associated density
    AgeDensity_N = [AgeDist_N(1, :); diff(AgeDist_N, 1, 1)];

    % Computes the function D_n^1(y) - see Equation (RP-34)
    Dy1 = Employ_Dist(p, F1, AgeDensity_N, Deltagrid, LPNwMLam);

    % Exogenous employment grid over which we want employment
    Dy_grid = exp(p.yvec)';

    % Creates a matrix for D_n^n(y) - see Equation (RP-31)
    Dy_n = zeros(p.ypoints, p.maxprods);

    % For each product we consider...
    for n = 1:p.maxprods

        % ... transform distribution according to Equation (RP-31)
        Dy_n(:,n) = interp1(n*Dy_grid, Dy1(:,n), Dy_grid);

    end

    % Remove NaN associated with off-grid interpolation
    Dy_n(isnan(Dy_n)) = 0;
    
    % Aggregate the distribution over the number of products - see Equation (RP-32)
    Employment = sum(Dy_n .* repmat(ProductDist, p.ypoints, 1), 2);

    % Make sure the resulting distribution adds up to 1 by correcting numerical imprecisions
    Employment = Employment/Employment(end);
    
    % Obtain the associated density for log(employment)
    Emp_density = [0; diff(Employment)] / (p.yvec(2)-p.yvec(1));

    % Transform density of log(employment) to density of employment
    Emp_density = Emp_density ./ Dy_grid;

    % ---------------------- Employment Shares ----------------------------

    % Determine the endpoints for the binned distribution
    emppoints = [4 9 19 49 99 249 499 999 2499 4999 9999];

    % Preallocate vector to store the share of employment in each bin
    empshare = zeros(1, length(emppoints)+1);

    % Denote previous employment point
    empold = 0;

    % For each endpoint we target...
    for i = 1:length(emppoints)
 
        % Integrand is simply the density-weighted employment
        Integrand = @(y) y.*interp1(Dy_grid, Emp_density, y, 'linear', 0);

        % Compute share of employment attributable to current bin
        empshare(i) = integral(Integrand, empold, emppoints(i)) / integral(Integrand, 0, Inf);

        % Update starting point for the next bin to be analyzed
        empold = emppoints(i);

    end

    % Computes the employment share of large firms - see Equation (RP-33)
    empshare(end) = integral(Integrand, empold, Inf) / integral(Integrand, 0, Inf);
    clear empold

    % Computes the employment share of large firms - see Equation (RP-33)
    ShareLarge = empshare(end);

    % ---------------------------------------------------------------------
    %                       5. Exit Rates by age                           %
    % ---------------------------------------------------------------------   

    % For an overview of this part of the code, please refer to Section 2.7
    % of the Replication Package

    %----------------------------------------------------------------------

    % Preallocate vector for the probability of exit given number of products
    exitrates_prod = zeros(1, p.maxprods);

    % Grid for the number of products that could be destroyed in a period
    Kprods = 1:p.maxprods;

    % For every possible number of products that a firm may have...
    for n = 1:p.maxprods

        % Vector that stores the terms to be summed over m
        Finalterm = zeros(1, p.maxprods + 1);

        % For every possible number of products could be created
        for m = 0:p.maxprods

            % Compute the term of the outer sum - see Equation (RP-40)
            outer = exp(-p.xstar*n)*(p.xstar*n)^m/factorial(m);

            % Compute the term of inner sum - see Equation (RP-40)
            inner = exp(-(p.tau+p.d0)*n)*((p.tau+p.d0)*n).^Kprods./factorial(Kprods);

            % Add the terms of the inner sum - sum starts at n+m - see Equation (RP-40)
            inner = sum(inner(n+m:end));

            % Multiply the outer term with the inner sum - see Equation (RP-40)
            Finalterm(m+1) = outer*inner;

        end

        % Stores probability of exit for a firm with n products - see Equation (RP-40)
        exitrates_prod(n) = sum(Finalterm);

    end

    % Transforms exit conditional on number of products into exit conditional on firm age - see Equation (RP-41)
    Exit_age = repmat(exitrates_prod', 1, TimePeriods);
    Exit_age = sum(Pmatrix .* Exit_age, 1);
    clear inner outer Finalterm Kprods m;

    % ---------------------------------------------------------------------
    %                                                                     %
    % ---------------------------------------------------------------------    



    fprintf('-------------------------- Counterfactual stage ended -----------------------\n')
    fprintf('\n\n')


end